Electronic States of Graphene Grain Boundaries 
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We introduce a model for amorphous grain boundaries in graphene, and find that stable structures 
can exist along the boundary that are responsible for local density of states enhancements both 
at zero and finite (~ 0.5 eV) energies. Such zero energy peaks in particular were identified in 
STS measurements [J. Cervenka, M. I. Katsnelson, and C. F. J. Flipse, Nature Physics 5, 840 
(2009)], but are not present in the simplest pentagon-heptagon dislocation array model [O. V. Yazyev 
and S. G. Louie, Physical Review B 81, 195420 (2010)]. We consider the low energy continuum 
theory of arrays of dislocations in graphene and show that it predicts localized zero energy states. 
Since the continuum theory is based on an idealized lattice scale physics it is a priori not literally 
applicable. However, we identify stable dislocation cores, different from the pentagon-heptagon 
pairs, that do carry zero energy states. These might be responsible for the enhanced magnetism 
seen experimentally at graphite grain boundaries. 



I. INTRODUCTION 

Grain boundaries and other extended defect structures 
in graphite have been studied by surface measurements 
techniques for quite some time .'iJ'^ This research actually 
reaches beyond the fundamental questions of mechanical 
material properties and crystalline ordering complexities. 
The study of defects on the surface layer of graphite are 
directly related to the influence of disorder on isolated 
graphene sheets, and thereby of direct relevance in the 
context of graphene's extraordinary properties and po- 
tential electronic applications. Grain boundaries have a 
special status, since they are the natural extended de- 
fects also in two-dimensional graphene, while they have 
a topological status since in terms of the lattice order 
they can be represented as an array of dislocations with 
Burgers vectors that do not canceLlS 

The STS studies of graphite have also revealed some 
clues about the connection of extended defects and the 
controversial ferromagnetic properties of metal-free car- 
bon.'Sl Earlier theoretical studies aiming at localized de- 
fects in graphene,l2EI do yield some insights into the elec- 
tronic states and magnetic properties of some types of 
graphene edges, cracks, and single atom defects. How- 
ever, theoretical studies of the extended defect structures 
themselves have been completely absent until recently.^ 

The recent STM and STS studies of the electronic 
properties of defect arrays in graphite^ have shown that 
the local density of states (LDOS ) has two types of char- 
acteristic features: either an enhancement at zero energy, 
or a pair of peaks at low energy below symmetrically 
distributed around the Fermi energy. A first-principles 
model of grain boundaries based on a periodic array 



of the simplest pentagon-heptagon dislocationa^ revealed 
the possibility of forming bands around zero energy when 
the dislocations are close to each other, accounting for the 
LDOS peaks at finite energies. 

Motivated by the STS measurement results, we aim 
at extending the theoretical knowledge of extended de- 
fect structures by analyzing the electronic structure of 
amorphous tilt grain boundaries in graphene, expecting 
our results to be directly applicable to measurements on 
the surface of graphite. Our approach is based on con- 
sidering the relaxed boundary of misaligned grains of 
graphene, and the results should be of direct relevance 
to the structures found along the grain boundaries as 
seen on the surface of graphite.'^^ We find that the dis- 
ordered structures formed at the relaxed boundary be- 
tween two differently oriented grains can have enhanced 
LDOS at zero, or at finite energies. These features result 
from narrow bands (localized states) that can form both 
near and away from zero energy. We also discuss grain 
boundary models derived from dislocation arrays, consid- 
ering dislocation cores that are different from the simplest 
pentagon-heptagon structure.'^ These can lead to LDOS 
enhancement at zero energy as seen in the STM measure- 
ments, that are not seen in the pentagon-hexagon model 
of Ref. l5i Finally, we do identify a special limit where 
the zero modes of the low energy continuum theory of 
dislocated graphene precisely agrees with tight-binding 
model results. Intriguingly, this theory predicts the ap- 
pearance of localized zero energy states in an array of 
well separated dislocations, in contrast to the results of 
the first principles calculations of Ref. 5 

This paper is organized as follows. In Section In] we 
review the LDOS of dislocations, considering two defects 
at different distances, as well as the isolated case. Next in 
Section [II B| we use the continuum theory of graphene to 
explain the density of states and predict the zero energy 
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peak in an array of dislocations. In Section III we present 
our study of the tight-binding model of relaxed tilt grain 
boundaries in graphene, with a variety of opening angles. 



We close with discussion and conclusions. 



II. DISLOCATIONS IN GRAPHENE AS BASE 
OF GRAIN BOUNDARY MODELS 

Dislocation models of grain boundaries rely on the fact 
that an array of dislocations with same Burgers vectors 
produces a boundary line between two crystal domains 
of different lattice orientationsl^^^ In this Section we 
make observations relevant to such models in graphene, 
inspired by the recent STM experiments.'^ 



A. Graphene dislocations in tight-binding 

Simple dislocation cores in graphene come in two 
shapes that we label "PH core" and "OCT core" (cf. 
Fig. IlTa),(d)), both of wh ich were shown to be stable 
lattice configurations .liiES Geometrically, the two possi- 
bilities arise because the Bravais lattice has two atoms 
(separated by A) in the unit-cell, so that there are two 
inequivalent mutual configurations of A and the Burg- 
ers vector b. The L DOS at the atoms forming the core 
has been considered,'^^^^ revealing a sharp peak at zero 
energy in case of the "OCT" core, due to the undercoor- 
dinated atom. Note that even if only tt— orbitals are con- 
sidered, the single excess atom in one sublattice carries an 
LDOS peak and the accompanying (locally unb alanced) 
magnetic moment in the presence of interactions! ^^ ^ ^^"^^ ^ 
Alternatively, the "OCT" core can be viewed as a piece 
of a zigzag graphene edge of minimal length of one atom 
embedded in the graphene bul k, le ading to same conclu- 
sions about its LDOS features .1212211211 

Ref. 5' considers only "PH" type cores as building 
blocks of grain boundaries, and such models have been 
proposed in earlier graphite STM measurements. GI^Sl 
However, the zigzag oriented grain boundary model of 
Ref. 6j as well as simple geometrical considerations we 
present above, both show that the arrays of "OCT" type 
dislocations should not be disregarded in real materials, 
even if they are more energetically costly than the "PH" 
type. 

The inclusion of the "OCT" dislocations can be im- 
portant for explaining the observed LDOS peaks at zero 
energy in the measurements of Ref. [BJ The set of grain 
boundaries considered there shows that such LDOS fea- 
tures are found only when the defect cores are well sep- 
arated (i.e. the grain boundary angle is small). One 
might assume that when the defects are closer to each 
other, the zero energy states hybridize and move to fi- 
nite energies. We have however found that the localized 
zero energy modes are robust even when the defects are 
brought next to each other, which would be the case in 
a grain boundary with maximal opening angle. 

Our analysis was done by considering the LDOS of 
defects set inside a 75x75 unit-cell sized graphene patch 
tight-binding model, with twisted periodic boundary con- 



ditions in both directions. The special boundary condi- 
tions enable the system under consideration to actually 
be a periodic, 10x10 sized arrangement with the graphene 
patches as unit-cells, thereby leading to a tenfold in- 
crease in linear system size and a correspondingly denser 
energy spectrum £"„ (with corresponding eigenfunctions 
ijjn): from which the LDOS p{i,E) at site i and energy 
E follows in a standard way: 



p{i,E) 



^Ei^"«i^i-^ 



E- En + ie 



(1) 



We applied a small broadening e « 20 mcV of levels into 
a Lorentzian shape, which is both expected to exist in the 
material and leads to smoothing of the finite size effects 
in the LDOS . We introduce the defects by inserting a 
line of extra atoms, thereby creating a defect — anti-defect 
pair at a maximum separation of half the graphene patch 
size. By adding an additional line of atoms, we can study 
the LDOS of two defect cores close to each other, isolated 
from their anti-defects. The Hamiltonian is of the single 
particle spin degenerate tight-binding graphene: 



^ = -E^ 



U V^i ^j 



H.C.), 



(2) 



<ij> 



with the hopping constant t — 2.7 eV. When choosing 
the nearest neighbor pairs in Eq. ([2]), we retain the topol- 
ogy of the honeycomb lattice, which is violated only at a 
single atom in the "OCT" dislocation case. The LDOS 
turns out to be robust to relaxation of bond lengths, so 
that the results for tij = t are representative. 

Our calculation shows that the LDOS at the disloca- 
tion cores is insensitive to the distance between the dis- 
locations, in particular the LDOS peak at zero energy 
in the "OCT" type core system stays pinned and does 
not hybridize when the defects are brought close to each 
other to minimal distance of few lattice constants. 

The results of the tight-binding model presented in this 
section show that the characteristic features of the dis- 
location LDOS , notably the zero energy peak of the 
"OCT" core fall-off with distance from the core as a 
power law (Fig.IlTe)). This is the expected behavior ac- 
cording to low energy continuum models of graphene (see 
Section II B I , and also argued for in the case of cracks in 



graphene in Ref. 7, 

The STS measurements of Ref. i6| achieving atomic res- 
olution, however show an exponential fall-off of LDOS 
features with the distance from the prominent defect cen- 
ters, even for defects far from each other. This discrep- 
ancy might be due to subtle shortcomings of substitut- 
ing a simplified single graphene sheet for the top layer of 
graphite; however, another explanation could be the pres- 
ence of stronger disorder. The fact that the single atom 
resolution along the grain boundary is lost in patches 
of several lattice constants across also indicates that the 
grain boundary might contain more disorder than an ar- 
ray of simple dislocations. This presents additional mo- 
tivation for our study of amorphous tilt grain boundaries 
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FIG. 1: The LDOS of graphene dislocations from tight-binding and continuum theory, (a) LDOS of representative 

atoms of the "PH" type dislocation core (inset, see Section II A). As discussed in section IIB, the weight is shifted 

due to the A-A bond (thick in the inset) , compared to the symmetric curves in (c) obtained by switching off the 

A-A bond that are consistent with continuum theory . (b) The influence of the A-A bond on an isolated ring: the 

lattice case (a) can be viewed as a broadened version, (c) "PH" core LDOS without the AA bond. Dislocation 
topology effects from continuum theory (Section IIB dashed black curve) are prominent. (Finite LDOS at i? = is 



a finite size effect.) (d) LDOS of representative atoms in "OCT" core type, (e) The height of LDOS peak in (d) 

(pentagon, green) falls of like a power-law with distance from defect core. This holds also for LDOS features in (a). 

(f) The continuum defect topology prediction of LDOS (in patch of radius 6 = 0.1, Section IIB I in a dislocation 

array, with zero modes. 



presented in Section III in place of the coherent ones 
studied in Refs. .1^5^26. 



B. Continuum model of dislocations 

It is interesting and fundamental to approach the de- 
scription of grain boundaries by considering an analytical 
model. In this Section, we describe the results of such a 
continuum model, finding conditional agreement with the 
tight-binding results. We then proceed to use the theory 
for describing the LDOS of an array of dislocations, and 
find a surprising prediction of localized modes at zero 
energy. Even if the continuum theory prediction fails in 
a more realistic model (as Ref. ,5] suggests), we find it a 
fundamental step in understanding the system. 

The continuum description of the topological effect of 
dislocations is based on the description of the defect as 
a translation by the Burgers vector b of the wavefunc- 
tion of the ideal crystal, upon encircling the defect core. 
The model is therefore akin to an Aharonov-Bohm (AB) 
effect, except that it does not break time reversal sym- 
metry. The details of this model are derived in Ref. \Tl\ 
and here we start from the Hamiltonian in the form of 
the standard graphene Dirac equation, coupled to a dis- 
location gauge field. 



where the dislocation gauge field A (in fixed gauge) pro- 
duces the correct pseudoflux of the translation holonomy 
§A-dx=iK- b)T3 = 2nd rg, e.g. A^ = i^^rg = ^rg, 
where r and (p are the standard polar coordinates. The 
Burgers vector is encoded in the dislocation pseudoflux 
d which has only three inequivalent values {0^ , — ^ } = 
{0, — |, — §}, opposite at the two Fermi points.f^We label 
a Fermi wavevector by K (and the other Fermi point is at 
—K), T matrices mix the two Dirac points, the a matrices 
act on the A/B sublattice, and we use the four compo- 
nent spinor vE'(r) ee {^k+a, *x+s, *i<:_s, -"^k.aV- 
An important property of the translation operator, and 
consequently the A gauge field, is that it does not mix the 
Fermi points, so that we can consider them separately. 
Therefore our model is based on a single valley Dirac 
equation in the AB field of flux d E {— |,— |}, 



m 



-ihvpa ■ (V 






(4) 



The other valley experiences the complementary flux 
-1 - d, i.e. fff = H^^~'^. We have chosen the val- 
ues of d such to conform to the practice of AB flux being 
the fractional flux part. 

To test this theory, we find that the LDOS in a patch 
of radius 5 covering the defect behaves as. 



piS,E) 



(5) 



Hdisi = -ihvFTo (g) a • (V - iA), 



(3) 



This actually agrees with the tight-binding model results 
in the limit where the bipartiteness of the honeycomb 



lattice is not broken by the defect, see Fig. fT|c). This 
is precisely the limit where we expect that the effects 
of the global topology in the hopping network become 
dominant. This condition can be realized in principle 
for both cores pending their 'chemistry'. In the "OCT" 
case the undercoordinated atom appears as an intruder, 
but otherwise the bipartitness and topology of the ideal 
lattice are preserved. In the "PH" core case, the A-A 
bond (inset of Fig. [I^) spoils the hopping bipartiteness 
when it supports a finite hopping. It interferes with 
the purely topological effect of the dislocation, and in- 
troduces asymmetric features in the LDOS (Fig. flja)) 
while the powerlaw behavior expected from the contin- 
uum limit is recovered when the bond is switched off 
(Fig. fl|c)). The origin of the asymmetric features is 
clearly identified by considering the LDOS of an isolated 
10 atom ring which is turned into a pentagon — heptagon 
structure by switching on the A-A bond (Fig. flVb)): the 
lattice results (Fig. fija)) can be viewed as the 'molecu- 
lar' states of Fig. ^(h) turning into broadened, resonant 
impurity bound states. 

We now outline the calculation leading to Ec^. (l5]), 
which is also fundamental for understanding the predic- 
tion for a dislocation array. The eigenfunctions of Eq. Q 
are found by separating the angle, and we find for energy 
E = ehvpX (e = ±1 and A > 0): 



*! 



M = E 



^tm^ 






(6) 



where the sign s = +, — labels two linearly 
independent solutions ip^ = (u^(r), w^(r))^ — 
(Js(m_i_d)(Ar), Js(m-d)(Ar'))'^, and Jq is the Bessel func- 
tion of order q (note that q ^ Z). The total angular mo- 
mentum in channel m is j = m — 1/2, and we see that 
the presence of dislocation shifts it j ^f j ~ d. Normal- 
izability allows exclusively V'm ^'^^ m > 0, and ipz.ior 
m < 0, and both for m = 0. The Hamiltonian Eq. Q is 
actually not self-adjoint, so that there is additional phys- 
ical input needed regarding the wavefunction boundary 
condition at the singular point at the defect. At this 
point the theory becomes sensitive to the 'UV, that is 
the microscopic details at the lattice cut-off. In the field 
theoretical derivation that follows this UV regularization 
is kept as featureless as possible. However, we find from 
the explicit tight binding description that the 'chemistry' 
of the core structure does matter. Even without the A-A 
bond, which spoils the global topology, the PH core can- 
not be represented by the continuum theory (Fig. flTc)). 
However, a key result of this paper is that the OCT dis- 
location core (Fig. [lVd),(e)) is compatible with the con- 
tinuum theory zero mode (Fig. [l|f)). 

The application of the standard theory of self-adjoint 
extensions (SAEp^EH prescribes that the coefficients of 
the linear combination N^tpQ + N^ipQ in channel m — 
determine the additional physical parameter x G [0, 27r) 
through 



The channel m = actually contains normalized spinors 
■0™ which have diverging components on the sublattice 
A/B respectively, and the ratio of these divergences is set 
by the particular SAE through the value of x- 

We can now evaluate the LDOS in a patch of radius 

p{5,E) ^2n rdrVV|V^„(Ar)|25(S~^A)^ (8) 

Jo , X ™ 



£,A m 



\^1. 



iV+/7V_=cot(x/2). 



(7) 



AdA / rdr {rXy^ S{E - hvpX) ^ (9) 


_j2g+2|^|2g+l^ (10) 

where in the second line we have used the Bessel function 
density of states J^x ~^ I XdX, and evaluated the small 
argument (i.e. Xr <C 1) expansion of the Bessel func- 
tions of order q. The leading contribution comes from 
the diverging components of ipo, where q = —1/3, —2/3 
(this holds for both Fermi points, and both dislocation 
classes) . 

The value oi q — —2/3 generates an unphysical diver- 
gence p ^ 1/|£'|^/^. We therefore have to choose the 
SAE which removes the offending part of ■00) and this 
turns out to be x = tt and x = Oi for d = —1/3 and 
d — —2/3, respectively. Note that at the second Fermi 
point, we have to switch the values, so that x = 0,7r for 
d — -1/3,-2/3. The surviving components in 0o with 
q = -1/3 yield the advertised patch LDOS of Eq. ^. 



C. Continuum model of dislocation arrays 

Once we know the details of the continuum description 
of a graphene dislocation derived in the previous subsec- 
tion, we can ask the question: what happens in an array 
of such defects? As we have shown, the SAE of the con- 
tinuum Hamiltonian of Eq. Q is fixed by the allowed 
values of X = 0, tt. It turns out that precisely these spe- 
cial values of x allow the Hamiltonian to have localized 
states at zero energy. This leads to a peak at zero energy 
which is absent from the gapless, cusp shaped LDOS of 
the finite energy wavefunctions in Eq. ([5]). 

It is well known, that Hamiltonians with singular po- 
tentials (e.g. AB flux,^^! Coulomb potential, ^"^ delta func- 
tion potentiaPH) , once they are made Hermitian through 
a SAE, can exhibit finite or zero energy bound states, 
even if the original Hamiltonian was scale-free. 

Our system is represented by two copies (two Fermi 
points) of a two-component, two dimensional spinor in 
the presence of a pseudomagnetic solenoid with flux d G 
{—1/3, —2/3}. The problem of a spinful two dimensional 
particle moving in an arbitrary magnetic field, both non- 
relativistic (Pauli) and relativistic (Dirac), has originally 
been considered by Aharonov and Casher,'^^ who found 
that the number of flux quanta give the number of zero- 
energy states of the particle. In Ref.f551 the Dirac particle 
in the presence of multiple AB solenoids is considered, so 
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FIG. 2: The LDOS and states of tight-binding amorphous tilt grain boundary: example of armchair type with 

medium opening angle, (a) The system with hoppings (of different strength in calculation) included. (b),(c) 

Zoom-in of the boundaries, including the wavefunctions at fc^ = at energy of the LDOS peaks: size of colored dots 

is the amplitude, orange (dark grey) and yellow (light grey) denote opposite sign. (d),(e) The LDOS of two grain 

boundaries, averaged within square patches as marked in (a). 



we can here directly use those results concerning the zero 
modes of the Dirac Hamiltonian of the form Eq. Q. 

Let us describe the relevant calculation, following 
Refs. I33I34I closely. The fact that x = 0, tt is the key 
ingredient: as we have seen in Section II B at these 
values the divergence of the wavefunction is allowed in 
only one of the spinor components. This means that 
the SAE imposed boundary condition on the wavefunc- 
tion does not mix the two components, i.e. sublat- 
tices. For zero energy, the eigenproblem of H^ also de- 
couples the sublattices. Going to complex coordinates 
z = X + iy, A = Ax -f iAy, and using the scalar potential 



^(z) = — ^" di log \z — Zi\ of the AB gauge potential of 
dislocations di at positions Zi (i.e. dz^{z) = A), we can 
solve for the two sublattices separately: d^* (e~*w) = 



and dz{e 



0. The Dirac equation now tells us that 



e~ u (e w) is an analytic (antianalytic) function outside 
the singular points Zi. For x = tt, m cannot have sin- 
gularities according to Eq. ([7|. Taking into account the 
behavior e~* ~ \A'^ ^ kl ~^ o°; with = —'Y^^di the 
total pseudoflux, it follows that u can be a polynomial 
of z of order at most {—</>} — 1, where {} is the lower 
integer part. There are {—</>} linearly independent such 
polynomials. In the case x = 0, w is not singular so that 



e*ti vanishes at the defects, t; is a polynomial in z* of 
degree {</>} — 1, with n zeros, and there are {4> — ri\ of 
them. 

Collecting the results, there are {|0 — n\} {{\(j)\}) zero 
modes in the case X = (x = "") for the single Fermi 
point system. Note that we assumed the same value of 
X for each dislocation di, so that the result holds only 
in the case of all dislocations having equivalent Burgers 
vectors, which is the case of a grain boundary! The two 
Fermi points contribute independently to the number of 
zero modes, and since x and di are reversed between 
them, we get a total of 



lattice onto the torus (i.e. the plane with periodic bound- 
ary conditions), but this is possible only when for every 
dislocation there is an anti-dislocation. In that case, the 
theorem predicts no zero modes, in accordance with our 
present calculation. 



III. ELECTRONIC TIGHT-BINDING MODEL 

OF RELAXED SYMMETRIC AMORPHOUS 

GRAIN BOUNDARIES IN GRAPHENE 

A. The Method 



D^{\2c^-n\} 



{|0|}, with0=-. 



(11) 



zero modes in graphene with an array of n disloca- 
tions having the same Burgers vector (of whichever non- 
trivial class d). This number takes the values D — 
2, 2, 2, 4, 4, 4, 6 . . ., starting at n = 4 and onwards. More 
precisely, we get D = 2{n/3}, so that D scales with the 
system size, i.e. Z? ^ |n in the thermodynamic limit. 

The zero energy modes are localized at the defects 
and have a power-law shape. To answer the question of 
whether they are observable, we look at how the LDOS 
in a patch of radius 5 scales in comparison to the LDOS 
contribution of the finite energy states Eq. (Is]) . Near the 
defect at z^, at one Fermi point the u spinor component 
(sublattice A) scales as \z — Zi\~'''{z — ZiY , where p < 
{—0} — 1 and the value of d is —1/3. This gives a contri- 
bution p{5) ~ j2p-2d+2^ rp-j^g same sublattice at the other 
Fermi point contributes through v ^ \z — Zi\^'^'^{z — ZiY, 
with t < {n — (j)} — 1, giving p{S) ^ (52t-2d^ ^^ ^-j^g ^^^gg 
of the opposite defect type, we get the same scaling, but 
on sublattice B. The leading contribution in the LDOS 
comes from the minimal values of p = and t — 0, giving 
one mode at the defect with the LDOS 



p(")(,5)^52/3_ 



(12) 



The scaling shows that the zero mode contribution p^^' 
is more favourable than the finite energy contribution 
p{S,E) ~ S^^^E^/^ at smaller (5, because the strongest 
zero modes are more localized than all the finite energy 
wavefunctions . 

Since in this Section we dealt with a Dirac particle, 
it is interesting to consider the number of zero modes 
through the Atiyah — Singer theorem: In graphene with 
disclinations this was already analyzed in Ref. [35] by us- 
ing the defect gauge field of a disclination. There it was 
shown that the number of zero modes is proportional to 
the Euler characteristic of the manifold. The key to the 
application of the theorem is that graphene with discli- 
nations can form compact manifolds, e.g. the fullerene 
molecule, so that the mapping of the lattice and hopping 
topology onto a compact continuum manifold is correct, 
leaving the low energy Dirac particle description valid. 
For the case of dislocations however there is no such pos- 
sibility. The exception is the map from the dislocated 



We consider two misoricnted graphene grains of same 
width, confined in a periodic box of width Ly and length 
Lx- The nominal box boundaries at y = 0,Ly are po- 
sitioned through the middle of the width of the "first" 
grain. The second grain is generated in the middle half 
of the box (from Ly/A and 3Lj,/4) and both grains are 
periodic with the box in the x direction. The bound- 
aries at a; = and x = L^ have twisted PBC so that 
the system is a periodic crystal of length N * L^ (we set 
N = 18), with momenta kx — 2tt/{NLx). The lattice of 
each grain is generated from its center, and terminated 
at the grain boundaries. The boundaries are symmet- 
ric, but in general the two grain boundaries do not have 
the same structure because the two grains have different 
centers of inversion symmetry. 

The allowed values for the grain's orientation 9i fol- 
low from the constraint of its periodicity with the box 
along the x direction. If L' = aiCi + bie2 is the vector in 
the basis of the graphene Bravais lattice which is to be 
wrapped along the x box direction, the constraint is 



cos ( 



6? 



h/2 
+ a,bi/2' 



as also explicated in Ref. [3S1 If we allow slight strain 
in the grain, the number of available orientations can 
be enlarged.*^ The grain opening angle 9 — 61—62 — 
2dx spans the entire [0°,30°] range (for both zigzag and 
armchair typ#I^. 

We relax the atoms in the system at zero tempera- 
ture using the molecular dynamics method, where the 
interatomic p otenti al for carbon is taken in the Tersoff — 
Brenner form.EZED xhe potential between atoms «,j at 
distance r,,- is 



VR{r)= ^ 
VA{r) - 



BijVAinj) 
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with /(r) the smoothing function 
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FIG. 3: The LDOS and states of tight-binding amorphous tilt grain boundary: example of zigzag type of medium 

opening angle, (a) The system with hoppings (of different strength in calculation) included. (b),(c) Zoom-in of the 

boundaries, including the wavefunction at fc^; = and E — 012 eV: size of colored dots is the amplitude, orange 

(dark grey) and yellow (light grey) denote opposite sign. (d),(e) The LDOS of two grain boundaries, averaged 

within square patches as marked in (a). 



The effect of bond angles is encoded in Bij — l/2{Bij 
Bji), with 



B^3 = 1 + E G(%fc)/(r,fc) 

where Oijk is the angle between the i— j and i — k bonds, 
and the function G is 



G{0) = «o 1 + 4 



dl d2 + (l+COs(0))2y 

The ground state of this non-spherically symmetric po- 
tential is the graphene honeycomb lattice, when the pa- 



rameters are chosen as in Ref. [351 D = Q eV, R = 0.139 
nm, /? = 21 nm-\ S = 1.22, 5 = 0.5, aq = 0.00020813, 
Co = 330, and do = 3.5. The smoothing cutoffs are cho- 
sen to include the nearest neighbor atoms, Ri = 0.17 nm 
and i?2 = 0-2 nm. 

When the lattice is formed, we consider a tight-binding 
model for electrons, of the form Eq. (pi). The hopping 
constants i^j are taken to fall-off exponentially, and fitted 
so that tij for the nearest neighbor distance |A| is i = 2.7 
eV, and for the next-nearest neighbor distance -\/3|A| it 
is t' = 0.1 eV, in accordance with accepted values for 
graphene. ^^ Finally, we extract the energy bands E{kx), 
the wavefunctions V-'b(*) and the LDOS p{i,E). 



B. Summary of Results 

We analyze in detail a number of grain boundaries 
of both zigzag and armchair type,'^'^ covering the en- 
tire range of opening angles by varying the box size: 
2.6a < Lx < 16.1a with a = 0.246 nm the graphene 
Bravais lattice constant. 

In summary, we find that the LDOS along the grain 
boundaries, averaged in square patches of size 4a and 
considered in the low energy regime of |-B| < 1 eV, shows 
three typical behaviors: 

(i) A peak at very small energy, \E\ < 0.05 eV 

(ii) Two peaks at nearly opposite energies, at around 
0.3 eV < \E\ < 0.5 eV 

(ill) Just one peak, at an energy 0.3 eV < \E\ < 0.5 eV. 

Focusing on case (i), we have determined that the 
lowest energy wavefunctions are sometimes localized on 
structures that resemble short zigzag edge segments (i.e. 
of length 2a) . This however occurs also in armchair type 
boundaries, but of course then the short zigzag segment 
is tilted away from the grain boundary line, the x axis. In 
some systems however, the zero energy peak is associated 
with overcoordinated atoms, having even five neighbours. 

We find that clear examples of case (ii) mostly appear 
at high opening angles (i.e. small L^), where the strong 
LDOS signal spans the entire grain boundary. There are 
also just a few energy bands with |_E| < 1 eV, so it is 
easy to identify that the LDOS peaks are due to van Hove 
singularities, in accordance with the findings of Ref. 'S'for 
large opening angles. 

The case (iii) we find is strongly correlated with car- 
bon atoms that were annealed into a position with four 
neighbors, meaning that four atoms are within the |A| 
distance, distributed roughly evenly around the central 
atom. 

Since the LDOS behavior of case (ii) has already been 
identified in Ref. IS] we illustrate the occurrence of cases 
(i) and (iii) through typical examples Figs. 2|3|4 



Finally, we note that typically there is one localized re- 
gion within the box that has atoms with high LDOS val- 
ues, i.e. one can say that there is one prominent "defect" 
within one L^. long unit-cell of the entire grain bound- 
ary. This means that the periodicity of grain boundary 
calculated from the opening angle corresponds to the pe- 
riodicity of prominent defect structures, even in our case 
of amorphous boundaries.*^ There are rare special cases 
where our box has an accidental symmetry so that the 
defect structures along the boundary repeat twice within 
the box length L^, effectively halving L^ and 9. 



Quite likely the grain boundaries that are formed spon- 
taneously in graphite, and that are best characterized 
experimentally, are of the relaxed amorphous kind as dis- 
cussed in Section III. Because of their disorderly struc- 
ture it is impossible to identify sharp and precise features 
in their electronic properties. Nevertheless, we do find 
that generically these support narrow bands at the grain 
boundary both close and away from the Fermi energy, 
of the kind seen in the tunneling experiments. A next 
question is what happens when the interaction between 
electrons is switched on. Due to the LDOS enhance- 
ment, we expect magnetic moments localized along the 
grain boundary, to be compared to the results of AFM 
scans of graphite in Ref. |6l This might provide a concrete 
model for (existence of) ferromagnetism found in defects 
on the graphite surface. 

We also analyzed in detail the electronic signature of 
ideal grain boundaries formed from arrays of dislocations. 
Starting from the perspective of continuum field the- 
ory revolving around the zero modes associated to Dirac 
fermions subjected to topological defects, we identified a 
potentiality of very elegant physics associated with grain 
boundaries. Combining dislocations in a grain bound- 
ary, we obtain the striking result that there are localized 
zero modes decaying as a power-law from the defect, and 
contributing to the observable LDOS . As we demon- 
strated, the relevancy of these field theoretical results 
are critically dependent on the details of the microscopic 
structure of the dislocation core. Murphy's law gets in 
the way with the most elementary and natural pentagon- 
hexagon dislocation core, possibly because this disrupts 
the topology underneath the continuum limit by spoiling 
the connectivity of the sublattices. However, the "OCT" 
dislocation core appears to be compatible with the con- 
tinuum theory, and we do find a zero mode structure and 
the correct power law behavior of the LDOS . 

We hope that our results will stimulate further exper- 
imental research. The challenge appears to find out how 
to control with great precision the microscopic structure 
of grain boundaries in graphene in the laboratory. In 
particular, it would be wonderful when it turns out to be 
possible to engineer grain boundaries formed from OCT 
dislocations, since this would form an opportunity to get 
a closer look at the profound beauty of the zero modes 
of Dirac fermions. 
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FIG. 4: The LDOS and states of tight-binding amorphous tilt grain boundary: example of armchair type of small 

opening angle, (a) The system with hoppings (of different strength in calculation) included. (b),(c) Zoom-in of the 

boundaries, including the wavefunctions at kx = at E of the peaks: size of colored dots is the amplitude, orange 

(dark grey) and yellow (light grey) denote opposite sign. (d),(e) The LDOS of two grain boundaries, averaged 

within square patches as marked in (a). 
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